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We present a generalization of the recently proposed variational cluster perturbation theory to 
extended Hubbard models at half filling with repulsive nearest neighbor interaction. The method 
takes into account short-range correlations correctly by the exact diagonalisation of clusters of finite 
size, whereas long-range order beyond the size of the clusters is treated on a mean-field level. For one 
dimension, we show that quantum Monte Carlo and density-matrix renormalization-group results 
can be reproduced with very good accuracy. Moreover we apply the method to the two-dimensional 
extended Hubbard model on a square lattice. In contrast to the one- dimensional case, a first 
order phase transition between spin density wave phase and charge density wave phase is found 
as function of the nearest-neighbor interaction at onsite interactions U > 3i. The single-particle 
spectral function is calculated for both the one-dimensional and the two-dimensional system. 
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I. INTRODUCTION 

In recent years an increasing number of theoretical and 
experimental studies in condensed matter physics have 
focused on the description and understanding of quasi 
one and two dimensional strongly correlated electronic 
systems. Several fascinating properties of these materials 
are due to the competition between different phases with 
long-range order. High-temperature superconductivity 
in cuprates is one of the most famous examples which is 
not yet understood in a satisfactory way. Realistic mod- 
els that are used in this context consist of a kinetic part 
which accounts for the electron motion and an interaction 
part which is of the same order of magnitude. The sim- 
plest model that can be constructed under these assump- 
tions is the tight binding Hubbard model. It consists of 
a kinetic energy part, where the electrons can only hop 
between nearest neighbor sites and the Coulomb interac- 
tion U which acts only locally on each site. Although this 
model was used with great success for the description of 
a wide class of materials, there are interesting physical 
questions which require an extension. The inclusion of 
the nearest-neighbor Coulomb interaction, for example, 
is necessary for the study of inhomogeneous phases, such 
as the charge-density wave (CDW). This leads to the so 
called extended Hubbard model (EHM). 

But knowing the appropriate model for the descrip- 
tion of a material is only the first step on the way to 
understanding the physics. Already for the simple Hub- 
bard model without non-local Coulomb interaction, an 
exact calculation of static and dynamic properties is pos- 
sible in very special cases only and one has to be con- 
tent with approximate methods in general. For the in- 
teresting case where the Coulomb interaction U is of the 
same order of magnitude as the bandwith W, the conven- 
tional perturbative approach must fail. This is expected 
for weak-coupling perturbation theory but also for the 
complementary approach with exact treatment of the in- 
teraction part and perturbative treatment of the kinetic 



1 9 

energypiifii* 

Numerical methods are more promising, such as quan- 
tum Monte Carlo (QMC)f^ exact diagonalisation (ED), 
and density-matrix renormalization group (DMRG)i^ 
They are able to give essentially exact results ~ at least 
for limited system sizes or (DMRG) for the one dimen- 
sional case. Another non-pcrturbative approach is the 
mean-field method and, in the context of the Hubbard 
model, the dynamical mean field theory (DMFT),^ in 
particular. While the DMFT directly works in the ther- 
modynamic limit of infinite system size, it must be re- 
garded as a strong approximation since spatial correla- 
tions are neglected altogether. Cluster generalizations of 
the DMFT include at least short-range correlations via 
the exact treatment of a small cluster instead of consider- 
ing a single impurity only. Both, a reciprocal-space (dy- 
namical cluster approximation, DCA^) and a real-space 
construction (cellular dynamical mean field theory, C- 
DMFTi2i2iia) have been suggested. 

Essentially the same idea is followed with the cluster 
perturbation theory (CPT)fiiii2ii^ which is a cluster ex- 
tension of the strong-coupling expansion for the Hubbard 
model: The lattice is divided into small clusters which are 
solved exactly while the hopping between adjacent clus- 
ters is treated perturbatively. The lowest order of the 
strong-coupling expansion in the inter-cluster hopping 
yields the CPT. Short range correlations on the scale of 
the cluster are taken into account exactly, for instance by 
the Lanczos technique at zero temperature, while correla- 
tions on a scale larger than the cluster size are neglected. 
The CPT is a systematic approach with respect to the 
cluster size, i.e. the method becomes exact in the limit 
Nc oo, where Nc is the number of sites within a clus- 
ter. It allows for the calculation of the single-electron 
Green's function at arbitrary values of the wave vector 
k. This is a considerable improvement compared to stan- 
dard Lanczos calculations for small clusters, where only 
a few k points are available. The CPT has been success- 
fully used to describe spectral properties of the high-Tc 
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niaterialsfiaii^ii& and has already been extended to finite 
temperatures 

Recently a new method has been proposed which ex- 
ploits a general variational principle for the self-energy 
of a system of interacting fermions. This self-energy- 
functional approach (SFA)i& approximates the self en- 
ergy of the original system in the thermodynamic limit 
by the self energy of an exactly solvable reference sys- 
tem with the same interaction part. The self energy is 
varied by varying the single-particle parameters of the 
reference system. Choosing the reference system to be a 
cluster of finite size yields a non-perturbative and consis- 
tent cluster approach. It has been showniSi that within 
this framework the CPT as well as the C-DMFT appear 
as special approaches depending on the number of ad- 
ditional uncorrelated ("bath") sites taken into account: 
The optimum number of bath sites is actually a free pa- 
rameter which can be determined from the general vari- 
ational principle. It has been pointed outlS, that at least 
for one-dimensional models a large cluster without bath 
sites has to be preferred. The use of a reference sys- 
tem without bath sites represents a generalized CPT in 
which the single-particle parameters of the finite cluster 
are optimized according to the variational principle. This 
"variational CPT" (V-CPT) has successfully been used 
in a recent study for the investigation of the symmetry- 
broken antifcrromagnetic phase of the two-dimensional 
Hubbard model. 

So far a consistent formulation of the (variational) 
cluster-perturbation approach could be achieved for lat- 
tice models with on-site interactions only. The reason 
for this restriction is that within the SPA the reference 
system must be chosen with the same interaction as the 
original model. As detailed in Ref.lTsl this ensures that 
functionals given by the skeleton-diagram expansion are 
the same for both, the original and the reference model. 
In case of the EHM the interaction couples the different 
sites of the lattice. Thus there is no reference system with 
the same interaction which consists of decoupled subsys- 
tems of finite size. The motivation of the present paper is 
therefore to extend the ideas of the CPT and V-CPT to 
the investigation of the EHM including nearest-neighbor 
Coulomb interaction. It is shown that a mean-field de- 
coupling of the inter-cluster nearest-neighbor interaction 
yields a systematic and reliable cluster approach. 

The paper is organized as follows: In Sec. we give a 
short description of the V-CPT method, Sec. lIIII shows 
how to decouple clusters in the case of the EHM. In 
Sec. lIVI and Sec. El we present results for one and two 
dimensions, respectively. The conclusions are given in 
SecEH 



II. VARIATIONAL CPT 

Let us consider a system of interacting fermions on a 
lattice with Hamiltonian H, in general consisting of a 
single-particle part and an interaction part Hi. The 



lattice is then divided into clusters, where it is of crucial 
importance for the derivation of the method that those 
clusters are connected by Hq only. The Hamiltonian can 
then be written as 

7? = [iff (R) + ffi(R)] + ^ i?^^R,R'), (1) 

where R denotes the individual clusters, ' (R) is the 
part of the single-particle term that acts only inside a 
single cluster, _ffi(R) is the interaction part inside the 
cluster, and the inter-cluster hopping is given by 

ij(^R,R')=ET^aT'4.^R',^.' (2) 

where the hopping matrix ^ is non-zero only for hop- 
ping processes across the cluster boundaries. The indices 
a and b are general quantum numbers within a cluster, 
e.g. position and spin index, and c\^^ ^ creates an electron 
with quantum number a in cluster R. 

The quantity of interest is the single particle Green's 

function GR,a,R',(,(t^) = {{c-R,^a^ c^^, ^ ■ Using transla- 
tional invariance at the level of the superlattice vector 
R, the Green's function becomes diagonal with respect 
to the wave vector Q from the reduced Brillouin zone 
corresponding to the superlattice. The resulting Green's 
function in reciprocal space is a matrix Gq(iu;) with el- 
ements Gci^a,b{'^) and a, & quantum numbers within a 
cluster. 

Within the CPT approximation this Green's function 
GQ(a;) can be expressed in terms of Green's functions 
of the decoupled clusters G'(cj), again matrices in the 
quantum numbers a and 6, and the intcr-cluster hopping 
T^j^ by the expression 

Gq(c^)= [G'(c.)-1-Tq]^' (3) 
with the Fourier-transformed inter-cluster hopping 

TQ,,, = l5:T,^,^'e^Q(^-^''. (4) 

R,R' 

For the details of the derivation of the CP T formulas we 
refer the interested reader to Refs. and references 

therein. We want to mention that one can transform 
Eq. ||2J) into a Dyson-like equation 

GQ(a;) = (Gg)(c.)-i-S(c.))-\ (5) 

where Gq''((jj) is the free Green's function of the infi- 
nite lattice, and S(ti^) is the cluster self energy. In other 
words CPT consists of approximating the self energy of 
the infinite system by the self energy of a cluster of finite 
size. Note that CPT is based on the exact evaluation 
of small clusters without any self-consistency procedure, 
and thus does not allow for the occurrence of symmetry- 
broken phases. This restriction is overcome with the V- 
CPT methodJi^ 
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The observation underlying V-CPT is that the per- 
turbation term need not necessarily be restricted to the 
inter-cluster hopping term but can be any single-particle 
operator. For this reason one has a certain amount of 
freedom for the partitioning of the single-particle part 
of Eq. 1^. Namely, the Hamiltonian Eq. is invariant 
under the transformation 



<)(R)^i/r(R)+o(R) 



r(i) 



(6) 



with an arbitrary intra-cluster single-particle operator 



(7) 



a,b 



which can for instance be a fictitious symmetry-breaking 
field, thus allowing for broken symmetry already on a fi- 
nite system instead of only in the thermodynamic limit. 
In the non-interacting case, where Eq. Q is exactfi2ii^ 
the result is independent of the transformation Eq. 
but in the interacting case the result depends on the par- 
ticular choice. However, this dependence is not a short- 
coming of the method but can rather be used for an op- 
timization procedure. 

The question of what choice for A = Aa^b will opti- 
mize the results can be answered by the SFA.^^ It pro- 
vides a way to exactly evaluate the grand potential 
as functional of the self-energy S by restricting the do- 
main of the functional to a certain subspace S of trial 
self-energies. This subspace consists of all S which are 
exact self-energies of the reference system for different 
A. The reference system H' must have the same inter- 
action part as H and must be exactly solvable for any 
A. Throughout the paper we use a cluster of finite size 
as reference system. 

The cluster self-energy can be parametrized as S = 
S(A). Within the SFA, the optimal value of A is deter- 
mined from the stationary point of the function— 

n{A) = n'{A) 

T y tr In —rj^^ — 



LT^ tr ln(-G'(A,iw„)), 



(8) 



where i7'(A) is the grand potential of the reference sys- 
tem. The frequency sum runs over discrete Matsub- 
ara frequencies icun, L is the number of clusters or Q 
points, respectively, T gives the temperature, and bold 
symbols denote matrices in the cluster indices a and b. 
Note that the fraction in the second line in Eq. ||SJ) is the 
CPT Green's function. The single-particle parameters 
A can include all single-particle parameters of the orig- 
inal Hamiltonian or only part of it, as well as additional 
terms, e.g. a fictitious staggered field. The more param- 
eters are considered for the optimization problem, the 
larger is the subspace S of trial self energies. The actual 



choice and number of parameters depends on the problem 
under consideration. For more details of the derivation 
of the method see Ref.S Note that the V-CPT method 
is superior to the CPT approach, since within the con- 
cept of V-CPT one gets the standard CPT formulas by 
setting the number of variational parameters A to zero. 

A necessary condition for the applicability of the 
method is that the clusters are coupled by single-particle 
operators only. At this point it is easy to see that a 
straightforward application of the method to the EHM 
where the clusters are also coupled by Coulomb interac- 
tions is not possible. However, we will show in Sec. lIIII 
how one can decouple the lattice into clusters appropri- 
ate for the application of CPT even in the case of the 
EHM. 



III. DECOUPLING THE CLUSTERS 

We start from the Hamiltonian of the extended Hub- 
bard model 

where i, j indicate the position in the lattice, and for con- 
venience we use a constant value Vi^ = V for all nearest- 
neighbor bonds. According to Eq. Q we decouple the 
lattice into clusters yielding 

R 

^ [ijf (R,R') + 4^R,R') 



R,R' 



(10) 



where the first row includes only terms of a single cluster 
and the second row couples different clusters. By com- 
paring the second row with the corresponding term in 
Eq. ^ one can see that the term causing problems in 
the case of the EHM is the interaction term 



h^\r,r') 



V 



y ^^R^nR'J■, 



(11) 



which is of two-particle type. The symbol [ij] indicates 
that the sum runs only over bonds connecting nearest 
neighbors in different clusters. For nearest-neighbor in- 
teractions this means that the indices in [ij] must belong 
to the cluster boundaries of two adjacent clusters. For 
the application of the method derived in SecElthe cou- 
pling term must be of single-particle type, which can be 
achieved by a mean-field decoupling of the interaction 
term Eq. (|ll|l . Hence we get 



M 

nR'j)- 

[y] 



(12) 
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Due to the translational invariance with respect to the su- 
per lattice vector R, the mean- field parameters (njn) and 
{riR'j) are independent of R and R' and will be denoted 
by Xi and Xj , respectively. With these abbreviations we 
get 

^v'mfII^' R') = 

R,R' 

= F ^ ^ [uRiXj + nw]Xi - XiXj] 

R,R' M 

= ^ [uRiXj + URjXi - X,Xj] 

= E^v.mf(R)- (13) 

R 

The double sum over R and R' reduces to a single sum, 
because for fixed values of R, i, and j only one term of 
the sum over R' contributes due to the fact that two-site 
interactions couple at most two different clusters. One 
has to be careful in order to avoid double counting of the 
bonds [ij]. For instance, for a one-dimensional cluster of 
length N, Eq. ifT^ reduces to 

[^riAat -I- urnXi - XiXn] , (14) 

R 

because the only decoupled bond connects sites 1 and N 
of different clusters. 

By this mean-field decoupling, two parameters Xi are 
introduced for each decoupled bond, e.g. Ai and Aat in 
one dimension, and in general all these parameters A.^ 
are independent of each other. But as we will see below, 
the number of mean-field parameters A^ can be strongly 
reduced in special cases. 

The decoupled interaction Eq. is of single-particle 
type and can be included in the intra-cluster hopping 
term Hq ' (R) , leading to a modified intra-cluster single- 
particle term 

Hq'^\'R, Xi) = Hq'^\'R) + i7y''j^^p(R, Xi), (15) 

where we explicitly denoted the dependence on the pa- 
rameters Xi. After mean- field decoupling we finally get 
the Hamiltonian 

HuFiXd = [^o'^R- AO + ^(R) + ^v^(R)" 

R 

+ ^i/(^R,R'), (16) 

R,R' 

for which the method described in Sec.^is applicable. 

From the decoupling of the clusters we have got ad- 
ditional parameters A^ which are external parameters to 
the Hamiltonian Eq. H16|l and have to be determined in 
a proper way. For this purpose we propose two different 
procedures: 

(i) One can get the parameters from a self-consistent 
calculation on an isolated cluster. That means that one 



starts with a certain guess for the Ai, which are the ex- 
pectation values of the electron densities on sites i. Then 
the ground-state wave function of an isolated cluster is 
calculated, giving new values for the A^. In this step 
open boundary conditions (obc) are used in order to be 
consistent with the obc necessary for the calculation of 
the cluster Green's function in Eqs. and ©. These 
new values Ai serve as parameters in the Hamiltonian 
for the next determination of the ground state, and the 
whole procedure is iterated until convergence of the Ai 
is achieved. This procedure may work quite well for the 
EHM in the case of a first order phase transition between 
a disordered and an ordered phase, because (due to an 
avoided level crossing) the transition point, i.e. the criti- 
cal Coulomb interaction Vc, is almost independent of the 
cluster size.^^ For second order phase transitions we ex- 
pect that this method will not give satisfactory results, 
because here we face a discrepancy between the parame- 
ters calculated on the isolated cluster and the parameters 
that would give the optimal result in the thermodynamic 
limit. 



(ii) The shortcoming in the case of second order phase 
transitions can be overcome in the following way: As we 
show in App. ^ the self-consistent calculation of mean- 
field parameters is equivalent to the minimization of the 
free energy F. Since the relation fl — F — fjiN holds at 
T = 0, this minimization can be done at the same time as 
the optimization of the single-particle parameters A in 
the SFA formalism, and we can use Eq. ||SJ) for the deter- 
mination of the parameters Ai, too. Note that all quan- 
tities in Eq. ((Sjl which depend on the single-particle pa- 
rameters A are dependent on the mean-field parameters 
Ai as well. To keep the calculations simple we consider 
only half-filled systems, where it is sufficient to use only 
two different values for the Ai, namely A^i = 1 — (5 and 
A_B = 1 + (5 on sublattices A and B, respectively. Under 
this assumption we have only one mean-field parameter 
5, and the grand potential is 57 = r2(A,i5). The general 
procedure is now, that for each value of 5 the stationary 
point with respect to A has to be found as required by 
the SFA formalism, yielding a function Vl = ^{5). By 
finding the minimum of this function one can determine 
the optimal value for 5. 



Conceptually, the latter method (ii) of determining the 
mean-field parameters is superior to the procedure (i) de- 
scribed first as it uses information on the Green's func- 
tion in the thermodynamic limit for the calculation of 5. 
However, one has to keep in mind that for each choice of 
5 the Green's function G'(cl') of the isolated cluster has 
to be calculated many times to evaluate Eq. © which 
is much more time consuming than the self-consistency 
procedure on the isolated cluster. 
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FIG. 1: Schematic phasediagram of the one-dimensional 
EHM taken from R.efs. l2 J2!Ti27i The thick line marks the first 
order phase transition, and the dashed line marks U — 2V . 



IV. ONE DIMENSION 



The Hamiltonian of the one-dimensional EHM is given 

by 



i.(7 i 

-I- y ^n^ni+i - /i^n^. (17) 



A. First order phase transition 

For a first test of our method we studied the one- 
dimensional EHM at [/ = 8, which is well above the 
multi-critical point. The phase transition is then of first 
order without any BOW phase between SDW and CDW 
phases. As reference system H' according to Sec^I we 
used decoupled clusters of different lengths consisting of 
Nc = 8, 10, and 12 sites, respectively. For the determi- 
nation of the mean-field parameter 5 we used the method 
(ii) described in Sec. lIIII where 5 is calculated from the 
minimum of the free energy of the system. For the SFA 
optimization of the single-particle parameters, we had to 
choose a set A of parameters which are varied in the 
optimization procedure. In order to minimize the num- 
ber of relevant parameters we used results of a recent 
study of the Hubbard model. There it has been shown 
that at [/ = 8 the variation of the hopping in the clus- 
ter yields only minor changes that can be neglected, and 
that open boundary conditions have to be used. More- 
over, for the one-dimensional Hubbard model it has been 
pointed out that the use of a fictitious staggered magnetic 
field as a variational parameter gives a stationary point 
of the grand potential only for vanishing field, yielding a 
paramagnetic phase without long-range magnetic order. 
Since it can be assumed that these results are also valid 
in a similar way for the EHM, we did not use the hop- 
ping in the cluster or a staggered magnetic field in the 
optimization procedure. Here we studied charge-ordering 
effects and therefore we used as variational parameter a 
staggered field coupled to the charge densities given by 
Eq. lO with 

Aa,6 = e^a,6e'Q^% (18) 



Throughout the paper we set t as the unit of energy. 
Although this model was studied intensively in recent 
years, the ground state phase diagram is still under 
some discussioni^^i^^i^'^1^^1^^1^^1^'^ We use this model as 
a testing ground for our method, because many results 
are available for comparison. The chemical potential is 
/i = C//2 -I- 2V due to particle-hole symmetry at half- 
filling. 

In one dimension at half filling, the EHM shows 
an involved phase diagram including spin density wave 
(SDW), charge density wave (CDW), and bond order 
wave (BOW) phases as shown schematically in Fig.^ 
There is good agreement on the existence of the latter 
phase, but its extension in the J7-y-plane has not yet 
been clarified in detailiSLS^ Moreover, a multi-critical 
point occurs at some critical value where the phase 
transition into the ordered CDW phase changes from 
a second order transition at lower values to a first or- 
der transition at higher values. Recent QMC studies^ 
gave Um ~ 4.7 for this point in good agreement with 
(^-ology2S*2i and bosonization results^^ whereas DMRG 
gives a somewhat lower value Um ~ 3.7i2& 



where Q = tt is the wave vector of staggered ordering 
and £ is the staggered- field strength. The grand poten- 
tial obtained in this way is shown in Fig.[21at two values 
of the inter-site Coulomb interaction. For comparison, 
calculations without optimization of the staggered field 
are shown as dashed lines in Fig.[21 As one can see, the 
optimization gives only minor changes to U,{5). The opti- 
mal staggered-field strengths in these calculations varied 
between Copt = 0.0 at 5 = 0.0 and Copt ~ 0.05 at i5 = 1.0 
at both values of V . 

From the shape of Vl{5) one can directly infer the order 
of the transition. If three minima occur at (5 = and 
5 =^ zbf^cDWi it is of first order, whereas it is of second 
order if ^{5) has only two minima at (5 = ±(5cdw f^nd 
a maximum at i5 = 0. As on can easily see in Fig.|21 
we have clear evidence for a first order phase transition 
at J7 = 8 with an SDW minimum a.t 5 = 0.0 and two 
degenerate CDW minima at i5 = ±i5cdw- At = 4.1 
the SDW phase is reahzed, ri(0) < ri((5cDw)) whereas at 
V = 4.2 we have fl{0) > ^{5cdw) and the CDW phase is 
the stable one. Thus we can state that the critical value 
Vc for the phase transition is located between V ~ 4.1 
and V = 4.2. 
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FIG. 2: Grand potential f2 as a function of the mean-field 
parameter S at U = 8 calculated on a cluster with Nc — 8 
sites as reference system. Upper panel: V — 4.1. Lower 
panel: V = 4.2. Solid lines: With optimization of a staggered 
field. Dashed lines: Without optimization of a staggered field. 




FIG. 3: Ground state energy Eq, kinetic energy -Bkin, and 
order parameter rriQ-Q^ of the one-dimensional EHM at (7 = 8 
after finite size scaling. Lines are guides to the eye only. 



For a more accurate determination of the phase bound- 
ary Vc, we have calculated the grand potential at several 
values of V and cluster sizes Nc = 8, 10, and 12. In addi- 
tion to the grand potential and the ground state energy 
Eq = fl + jiNe with Ne the number of electrons in the 
system, we calculated the order parameter 

mcDW = 1^ E - ^'"^^^ ' (19) 

where Q = tt, iVc is the number of cluster sites, and 
the kinetic energy i?kin. Both properties can be ex- 
tracted from the spectral function AR_r,R',r',o-(w) — 

(-l/7r)Im((cj^^^^^;c^,_^,^^})S''*^ The kinetic energy is 
obtained after Fourier-transformatiorkiSiiiS£ 

R,R' r,r' 

(20) 

by 

2 f" 

Ekin = jY. dt^e(k)A(k,w), (21) 

with e(k) the dispersion of the non-interacting system. 
Within our approach it is necessary to use the Lehmann 
representation for the cluster Green's function with small 
but finite Lorentzian broadening a. Whereas the grand 
potential Eq. ((SJ shows only minor dependence on this 
broadening, the dependence of the order parameter and 
the kinetic energy is considerably larger and one has to 
do an extrapolation to tr = Oi^ Although the formal- 
ism applies to the thermodynamic limit, results show a 



finite size dependence due to the finite size of the clus- 
ters serving as reference system. We found that the order 
parameter exhibits the strongest finite-size effects, which 
were of the order m-cDW n =io/"^CDw n =12 ~ 1-02 at all 
values of V. Finite-size scaling to A^c = 00 is easily done 
and the results for the ground state energy extracted from 
the minimum of the grand potential, the kinetic energy 
and the order parameter are shown in Fig.|31 Our results 
should be compared to Fig. 10 of Ref.f2ll which shows ex- 
cellent quantitative agreement with a deviation of less 
than 2 % for the calculated quantities at all values of V. 
From our calculations we get Vc = 4.140(5), again in 
agreement with the previous studies Refs. l2 1112611 

In order to provide a complete picture of the method 
we also performed calculations with mean-field param- 
eters obtained by a self-consistent procedure on an iso- 
lated cluster, see method (i) in Sec. lIIII For instance 
for Ac = 12 and V ^ 4.1 one finds self-consistent solu- 
tions for S — and for Ssc — 0.832, which differs only 
slightly from the value extracted from the grand poten- 
tial, (5cDW = 0.822. For this reason the calculation of 
the ground-state energy, kinetic energy, and order pa- 
rameter using 6sc instead of Scow gives practically the 
same results as in Fig.O In the present case it is there- 
fore sufficient to calculate the mean-field parameter from 
an isolated cluster which is much faster than finding the 
minimum of the grand potential. 

Whereas the properties we have shown so far are well 
known for the one-dimensional EHM, we additionally cal- 
culated for the first time the spectral function by Eq. H20() 
for arbitrary wave vector k. In Fig.^results are shown at 
[7 = 8 and selected values of V with a reference system 
consisting of A'c = 12 cluster sites, and the mean-field 
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FIG. 4: Density plot of the spectral function A(k, to) of the one- dimensional EHM at (7 = 8, calculated on a cluster of size 
Nc = 12 with Lorentzian broadening a — 0.1. Darker regions represent larger spectral weight. Coulomb interaction V as 
indicated in the plots. White lines are fits to a Hartree-Fock SDW/CDW dispersion (see text). 



parameter S calculated self-consistently by method (i), 
see Sec. mil We want to mention that the 'striped' struc- 
ture, particularly visible in the regions marked by 'C in 
Fig. 21 occurs because the decoupling into clusters breaks 
the translational invariance of the system. 

The spectral function a,t V — 2.0 is very similar to 
the spectral function of the Hubbard model {V = 0)iSii^ 
with splitting of the low-energy band into a spinon and 
an holon band, which are marked in Fig.0] by 'A' and 
'B', respectively. This similarity could have already 
been expected based on the full Hartree-Fock solution — 
decoupling of all interaction terms in the Hamiltonian — 
where one has no dependence on V at all in the SDW 
phase. But this simple picture holds only away from the 
transition point Vc as can be seen in Fig.^jin the plot at 

V = 4.0. At this point, in the vicinity of the phase tran- 
sition Vc = 4.14, the gap is considerably smaller than at 

V = 2.0, a clear deviation from the Hartree-Fock predic- 
tion. This indicates that charge fluctuations become very 
important in this regime, which are completely neglected 
by the Hartree-Fock approximation, but are taken into 
account on the length scale of the cluster in our approach. 



TABLE I: Fitted values for the hopping matrix element tat, 
gap Afit, and gap Ahf of the full Hartree-Fock approximation 
at (7 = 8. 





ifit 


Aflt 


Ahf 


V = 0.0 


1.93 


2.24 


3.75 


V = 2.0 


2.11 


2.20 


3.75 


V = 4.0 


2.62 


1.29 


3.75 


V = 4.5 


1.86 


3.35 


4.80 


V = 6.0 


1.86 


7.29 


7.88 



But although we found this deviation, one can still see 
residuals of the splitting of the low-energy band, a sig- 
nature for spin-charge separation. For this reason we 
infer that spin-charge separation is present up to the 
transition point. The white lines in Fig.^ correspond 
to fits of the holon branch to a Hartree-Fock dispersion 
E{k) = ±a/A2 -He(k)2. The fitted values for the hop- 
ping matrix element tat and the gap Afit are denoted in 
Tab.Q] where we included the values at V = for com- 
pleteness. One finds that the gap Afit is almost constant 
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from F = to y = 2 and, as mentioned above, consid- 
erably decreases near the the phase transition {V = 4). 
The hopping matrix element tfn shows the opposite be- 
havior and increases when approaching the transition 
point from below. This is due to the fact that in the 
vicinity of Vc, doubly-occupied and singly-occupied sites 
become close in energy, which enhances the movement 
of the electrons. The actual value of the matrix element 
tfit is very large compared to the original value t — 1 in 
the Hamiltonian. A fit to the spinon band would give 
a smaller value closer to t = 1, but whereas fitting to 
the holon band is consistent over the whole range of mo- 
mentum vectors k, the spinon band is only present for 
k < 7r/2 for uj ~ fi < (and k > 7r/2 for — > 0, 
respectively). 

The spectral function in the CDW phase shows a qual- 
itatively different behavior. At V = 4.5 we found a gap 
considerably larger than in the SDW phase, and this gap 
increases very fast with increasing V, as can be seen in 
the plot a,tV — 6. Moreover, no evidence for spin-charge 
separation can be seen in the spectral functions. By com- 
paring the fitted value A^t with the Hartree-Fock solu- 
tion Ahf, one can see that the agreement at V = 4.5 
is better than aX V = A, and that it becomes still bet- 
ter with increasing V. For this reason we conclude that 
charge fluctuations which are neglected in the Hartree- 
Fock approximation play a minor role in the CDW phase. 



B. Second order phase transition 

So far all calculations were done at U = 8, where the 
system shows a first order phase transition. In the fol- 
lowing, we study the EHM at U = 3, where the model 
exhibits a second order transition into the charge ordered 
CDW phaseii2iiS& In this paper we do not consider the 
BOW, since it has been argued that the SDW-BOW tran- 
sition is of Kosterlitz-Thouless type.-* For an analysis of 
this type of transition the available cluster sizes are far 
too small and do not allow a clear distinction between 
SDW and BOW phase. 

We calculate the grand potential fl{S) in the same way 
as in Sec. lIV Al in order to determine S. The result of a 
calculation on a cluster consisting of iVc = 8 sites is shown 
in Fig.|31 One can easily see a striking difference between 
the grand potential atU — 8, Fig. 12 and atU = 3. In the 
latter case there is only a single minimum. It is located at 
6 — tor V < Vc- With increasing V the curve for il{S) 
becomes flatter in the region around 6 — and finally 
two degenerate CDW minima occur at S = ±(5cdw for 
V > Vc- Note that here S changes continuously when 
crossing Vc, whereas it shows a discontinuity in the case 
of a first order phase transition. 

We find that now it is indeed important to use a stag- 
gered field, Eq. (|18l) . as a variational SFA parameter. In 
Fig.|Sl results are shown with such an optimization (solid 
lines) and without (dashed lines). Whereas at = 1.6 
both calculations show only the SDW minimum at 6 — 0, 




0.8 1 



FIG. 5: Grand potential fl as function of the mean-field pa- 
rameter 5 at i7 = 3 calculated on a cluster with Nc = 8 sites 
as reference system. Upper panel: V = 1.6. Lower panel; 
V = 1.7. Solid lines: With optimization of a staggered field. 
Dashed lines: Without optimization of a staggered field. The 
arrow marks the CDW minimum at V = 1.7. 
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FIG. 6: Ground state energy i5o, kinetic energy -Ekin, and 
order parameter mQYmi of the one-dimensional EHM at [/ = 3 
for cluster sizes Nc — 8 (dotted), Nc — 10 (dashed), and 
Nc = 12 (solid line). 



they differ at F = 1.7 where the system should already be 
in the charge-ordered phasei2ii22iS^iS& Without optimiza- 
tion of the staggered field, we would still find the SDW 
minimum at 5 = 0, but with optimization the minimum 
shows up for a finite value of ^ = ±0.31 characteristic for 
the CDW phase. 

For the determination of the critical value Vc, we cal- 
culated the ground state energy Eq, kinetic energy i^kin, 
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FIG. 7: Density plot of the spectral function A(]i.,uj) of the 
one- dimensional EHM at [/ = 3 calculated on a cluster of 
size Nc = 12 with Lorentzian broadening a — 0.1. Darker 
regions represent larger spectral weight. From top to bottom: 
V = 1.0, 2.0, 3.0. White lines are fits to a Hartree-Fock 
SDW/CDW dispersion (see text). 



and the order parameter meow at several values of V 
shown in Fig.|Hl We performed no finite-size scaling like 
in Sec. lIV Al since we found that here the cluster sizes are 
too small for a systematic scaling. Different from Fig.O 
E\iin, "T-CDW, and the slope of the ground state energy are 
continuous across the transition point as required for a 



TABLE II: Fitted values for the hopping matrix element tat, 
gap Afit, and Ahf within the Hartree-Fock approximation at 
U = 3. 







Aflt 


Ahf 


V = 


1.38 


0.29 


0.93 


V = 1 


1.59 


0.29 


0.93 


V = 2 


1.40 


1.13 


2.12 


V = 3 


1.59 


3.71 


4.28 



second order transition. From the kinetic energy and the 
order parameter calculated on a cluster of size Nc = 12, 
we extract a critical value of Vc = 1.665(5), which is 
in good agreement with the critical value Vc ~ 1.65 ob- 
tained by QMC^^ and diagonalization methods, ^'^■^'^ and 
with Vc = 1.64(1) from DMRG calculations.^ The slight 
difference is likely due to remaining finite-size effects. 
Moreover we made use of a single variational parame- 
ter only, namely the staggered field Eq. H18|) . and it can 
be expected that including more single-particle parame- 
ters in the SFA optimization procedure would give even 
more accurate results. 

We would like to point out that in the present case 
of a second order phase transition, the most accurate 
way of calculating the mean-field parameter 6 is to find 
the minimum in the grand potential including SFA opti- 
mization of single-particle parameters. Calculations on a 
cluster of size Nc = 12 showed that without optimization 
the critical value would be Vc = 1.685(5). Compared 
to Vc = 1.665(5) this is further away from the values 
obtained by other methods as given above. Calculations 
with S obtained self-consistently on an isolated cluster are 
insufficient. In this case one would get Vc — 1.735(5) for 
the TVc = 12 cluster. This means that for a second order 
phase transition S should be determined by minimizing 
the grand potential, whereas for first order transitions 
the self consistent determination was sufficient. 

The spectral function A(k,a;) a,tV= 1.0, 2.0, and 3.0, 
which has not been calculated previously, is depicted in 
Fig.[71 Wc found that the spectral function at y = 1.0 
shows only minor differences to the spectral function of 
the Hubbard model {V = 0). The white lines in Fig.[7| 
are fits to a Hartree-Fock SDW/CDW dispersion. The 
parameters tat and Agt can be read off from Tab.HTl In 
the SDW phase at V = and F = 1.0, the gap Afit is 
constant. Similar to the case U = 8 the agreement be- 
tween Afit and Ahf is better in the CDW phase than 
in the SDW phase. The hopping parameter tat increases 
when approaching the phase transition from below, simi- 
lar to Tab.lU but the fitted values for ifit are considerably 
smaller than in the case U = 8. 



V. TWO DIMENSIONS 

The two-dimensional Hubbard model is one of the 
most intensely discussed models for strongly-correlated 
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FIG. 8: Possible tilings of the two dimensional square lattice 
into clusters that allow for staggered ordering: Nc = S (bot- 
tom right), Nc = 10 (bottom left), super cluster with Nc = 48 
(top). 

electron systems, especially in the context of high- 
temperature superconductivity. But different from the 
one-dimensional case, where many sophisticated meth- 
ods have been used to investigate the extended Hubbard 
model as described in Sec. II VI only few studies have been 
done for the two-dimensional EHM. One reason for this is 
that many modern methods such as DMRG or fermionic 
loop-update QMC are difficult to apply to more than 
one spatial dimension. However, within our present ap- 
proach, the extension to two dimensions is straightfor- 
ward. 

The two-dimensional EHM is defined by the Hamilto- 
nian 

+ V^n^nj - ^^ni, (22) 

where {ij) connects nearest neighbors and the chemical 
potential is ^i = U/2 + AV at half filling. Early QMC 
studies^S showed that this model has a SDW-CDW tran- 
sition similar to the one-dimensional case with transition 
point Vc ~ C//4. But due to numerical difficulties it was 
impossible to determine the exact position and the order 
of the phase transition. Calculations within the Hartree- 
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FIG. 9: Grand potential calculated on a cluster of size Nc — 8 
at f/ = 8, 1/ = 2.1 (upper panel) and [/ = 3, F = 0.76 (lower 
panel) . 

Fock approximation^^ showed two stable phases for the 
Hamiltonian Eq. ^ at half fiUing, the SDW and CDW 
phase, separated by a phase boundary ai Vc — U /4. 

For the application of the method presented in Sec.lTTl 
the two-dimensional square lattice has to be decoupled 
into clusters of finite size. Three possible tilings with 
different numbers of cluster sites Nc are shown in Fig.|Sl 
Some care has to be taken concerning the staggered or- 
dering. Whereas for clusters with Nc = 8 and Nc — 10 
shown in Fig.|Sl the staggered ordering indicated by open 
and full circles is commensurate over the cluster bound- 
aries, a straightforward decoupling into clusters of size 
Nc — 12 is not possible. As one can easily see, a super 
cluster with Nc — 48 consisting of four iVc = 12 clus- 
ters has to be constructed in order to take into account 
the staggered ordering correctly. The Green's function 
of the super cluster can be calculated by switching off 
the hopping processes that connect the single Nc = 12 
clusters, in other words on bonds across the dotted lines 
in Fig.|Sl This gives a block-diagonal Hamiltonian which 
can be treated by the Lanczos algorithm. The switched 
off hopping processes are then incorporated again per- 
turbatively, that means by including the corresponding 
hopping terms in the matrix T^^^ in Eq. . Note that 
here the vectors R and R' denote the super clusters and 
not the single Nc — 12 clusters. Of course there are 
many other possible tilings like the 4x3 cluster used 
in Refs. 12, 1^.31, but also in that case a super cluster of 
Nc — 24 has to be used. 

We start the analysis of the two-dimensional EHM with 
the determination of the order of the phase transition. 
For this purpose we use the Nc = 8 cluster shown in Fig.|S| 
and calculate the grand potential Q{S) in the vicinity of 
the transition point at t/ = 8.0 and U = 3.0 as described 
in method(ii) in Sec. lHIl Here we did not use a stag- 
gered field as variational parameter, because it does not 
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change the quahtative shape of ^l{5) (see Figs. |21 and [SJ 
and is therefore not necessary for the determination of 
the order of the transition. The resuh of this calcula- 
tion is shown in Fig.|51 At both values of U we found 
three minima, located at S — and 6 = ±(5cdw- This 
indicates a first order phase transition, different from the 
one-dimensional EHM, where at [/ = 3.0 the transition 
is of second order. We checked that this different behav- 
ior is not likely to be a finite size effect due to the small 
linear dimension of the two-dimensional iV^ = 8 cluster 
by calculating 11 (5) for the one-dimensional model with 
A^c = 4 which still shows clear evidence of a second-order 
phase transition at C/ = 3. 

The fact that the system shows first order transitions 
at both I] ~ 8.0 and U = 3.0 simplifies the subsequent 
calculations. As discussed in the previous section, one 
gets good results in the case of a first order transition by 
using a mean-field parameter S determined self consis- 
tently on an isolated cluster, as described in method (i) 
in Sec. lIIII This procedure is much faster than the calcu- 
lation of the grand potential for many values of 6, which 
makes it possible to use the Nc = 48 super cluster shown 
in Fig.|Hl We want to mention at this point that the cal- 
culation of the grand potential for the two-dimensional 
system is much more time consuming than for one dimen- 
sion because of the larger number of Q points required 
in Eq. For one dimension L w 40 is sufficient for 
convergence, whereas L « 500 is necessary for two di- 
mensions. Nevertheless it is of crucial importance to use 
a cluster as large as possible, because the ratio of bonds 
treated exactly to mean-field decoupled bonds increases 
with increasing cluster size, especially pronounced for the 
two-dimensional square lattice. After having determined 
the mean-field parameter 5 for the CDW phase self con- 
sistently, we also performed an SFA optimization of a 
staggered field, Eq. ((TH|l . 

A few more words have to be said about calculations 
in the SDW phase ((5 = 0). Recent studies^^ of the pure 
Hubbard model revealed that it is important to take into 
account the long-range magnetic order for the accurate 
description of salient features of the system. This can 
be achieved by using a staggered magnetic field as varia- 
tional parameter, given by 

Aa^fa = Ma,fcZae^^^% (23) 

with Za- = ±1 for spin projection a —1,1, and h the 
strength of the field. Additionally it was argued that 
due to the connection of the hopping parameter t and the 
magnetic exchange constant J, results could be further 
improved by letting the hopping in the clusters be of 
strength t' and optimizing the staggered magnetic field 
and t' simultaneously. Therefore we use 

Aa.b = M^.bZ^e^Q^" - r5(,b>' , (24) 

where the symbol 5(^ab)' is equal to one for nearest- 
neighbor bonds inside the cluster and zero otherwise. 
The field strength h and t = t' — t are the variational 
parameters in the optimization procedure. 



U = 8.0 U = 3.0 
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FIG. 10: Ground state energy Eg, kinetic energy Skin, and 
order parameter meow of the 2D EHM at f/ = 8.0 (left) and 
U — 3.0 (right). Calculations were done on a A^c = 48 super 
cluster. 



To sum up, the following steps are performed in the 
analysis using the Nc = 48 super cluster: (i) First we 
determine the mean-field parameter (5cdw in the CDW 
phase self-consistently on an isolated cluster and (ii) use 
a staggered field Eq. 1)181) for an SFA optimization pro- 
cedure, (iii) In the SDW phase {6 = 0) the staggered 
magnetic field Eq. (|23|) and the intra-cluster hopping t' 
are optimized simultaneously, (iv) After determination 
of the SFA variational parameters we calculate the quan- 
tities we are interested in. 

The results for the ground state energy, kinetic en- 
ergy, and order parameter are shown in Fig. 1101 At both 
U = 8.0 and U = 3.0, the behavior of a first order tran- 
sition can be seen, where the change in the slope of Eq 
is much stronger at U ~ 8.0 than at U — 3.0. This 
change at U — 8.0 is even more pronounced than for the 
one-dimensional model at U = 8.0. From Fig. llOl we can 
extract the critical value Vc of the phase transition by 
fitting Eq to a straight line in the vicinity of the tran- 
sition point, and for the Nc = 48 super cluster we find 
Vc = 2.023(1) atU = 8.0, and K = 0.770(3) at U = 3.0. 
These values of Vc are much closer to the Hartree-Fock 
result Vc = U/A than for one dimension. Within our 
approach we cannot clarify whether this is an intrinsic 
feature of the two-dimensional model or it is an artifact 
of the approximation due to the larger number of mean- 
field decoupled bonds. 

The SFA variational parameters in the SDW phase 
near the phase transition point are found to be almost in- 
dependent of the interaction V. At U = 8, the optimiza- 
tion resulted in i' « 1.1 for the intra-cluster hopping and 
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/i » 0.14 for the staggered magnetic field. The optimiza- 
tion of just one single parameter leads to t' ~ 1.03|;i=o 
and h « 0.12|t/=i, and the value of also differs signifi- 
cantly from the value obtained by the simultaneous opti- 
mization of t' and h. This means that due to the strong 
connection between the magnetic ordering and the hop- 
ping matrix element it is important to optimize t' and h 
simultaneously in order to get the best approximation for 
the physics in the thermodynamic limit. In the charge 
ordered phase the dependence of the variational parame- 
ter, Eq. (|18|l . on the interaction V is larger with e = 0.08 
at V = 2.01 and e = 0.22 at V = 2.1. A similar be- 
havior can be found at U = 3: In the SDW phase the 
variational parameters t' « 1.61 and h « 0.15 are almost 
independent of V. In the CDW phase we get e ~ —0.03 
a.tV = 0.76 and e = -0.18 at V = 0.84. 

Whereas the application of the magnetic staggered 
field exhibits the symmetry h —h, this is not the case 
for the staggered field Eq. H18|l . because the symmetry is 
already broken by the mean-field decoupling. We found 
no stationary point of il for finite h in the CDW phases. 

The spectral function at f/ = 8 in the SDW phase 
{V = 1.0) and in the CDW phase {V = 3.0) is shown in 
Fig. ^2 We found that the spectral function at y = 1.0 
is very similar to the spectral function of the Hubbard 
model {V = 0)^ One can see that the spectrum mainly 
consists of four features, two high-energy Hubbard bands 
and two low-energy quasi-particle bands, separated by a 
gap in the spectrum. The dispersion of these low-energy 
excitations in the SDW phase differs significantly from 
the Hartree-Fock shape shown as white lines in the upper 
panel of Fig. llll which does not account for the splitting 
into coherent low-energy bands and high-energy Hubbard 
bands. The fit parameters were ifit — 1-34 and Afit — 
2.51. The width of the coherent bands \uj{X) — oj(r)\ « 
1.25 is rather set by the magnetic exchange J, consistent 
with QMC calculations atV = oMt^ The black lines are 
fits to E)^ = ± [-A + J/2{ cos fca; -t- cos fcy)^] which ac- 
counts better for the dispersion of the low-energy bands 
than the Hartree-Fock dispersioni^ The fit parameters 
were A^t = 2.69 and Jfit = —0.63, which is in good 
agreement with the second-order perturbation theory re- 
sult J = -4tV(C/ -V)^ -0.57. In the CDW phase the 
white lines correspond to Hartree-Fock dispersions with 
fit parameters Agt = 7.69 and tat — 1.16, and different 
from the SDW phase they agree well with the excitations 
of A(k,w). 

Fig.ll2l displays the spectral function at [/ = 3 at in- 
teractions = 0.5 and V = 1.0, respectively. The white 
lines again correspond to Hartree-Fock dispersions with 
fit parameters tat = 1.06, A^t = 0.64 at V = 0.5, and 
tat = 1-07, Afit = 1.82 at y = 1.0, respectively. As in 
the case U = 8 the dispersion of the coherent low-energy 
bands in the SDW phase differs from the Hartree-Fock 
prediction, but in this case the deviation is much smaller. 
We did not find an accurate functional form in order to 
fit the low-energy excitations, but nevertheless we can ex- 
tract the value of J from the band with of the coherent 
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FIG. 11: Density plot of the spectral function yl(k,aj) of the 
two-dimensional EHM at U — S calculated on a A'^c = 48 su- 
per cluster with broadening a — 0.1. Darker regions represent 
larger spectral weight. Top: V — 1.0. Bottom: V — 3.0. 
White lines are fits to Hartree-Fock dispersions. For the 
meaning of the black lines a.t V = 1.0 see text. 



bands yielding J = -(l/2)|w(X) - w(r)| « -1.57. This 
value is again in good agreement with the perturbation 
theory result J = -W^/{U - V) = -1.6. 

We would like to mention that our results at U ^ 3 in 
the SDW phase are qualitatively different from QMC re- 
sults at U = 3, V = 0, and inverse temperature p = 3tr^ 
where the spectral function shows metallic behavior with 
no gap around the Fermi energy. This difference may 
be due to temperature effects or due to poor resolution 
of the Maximum-Entropy inversion of QMC correlation 
functions. 

At both U — 8 and U = 3, one can easily see that 
agreement of the Hartree-Fock dispersions with the low- 
energy excitations of v4(k, oj) is better in the CDW phase 
than in the SDW phase. In addition the gap Ahf cal- 
culated within the Hartree-Fock approximation is much 
closer to the fitted gap A^t in the CDW phase (e.g. 
Ahf = 7.76, Afit = 7.69 at C/ = 8, V = 3) than in 
the SDW phase (e.g. Ahf = 3.57, Agt = 2.51 at U = 8, 
V = 0). Therefore we conclude that in the CDW phase 
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FIG. 12: Same as Fig-.E) but at [/ = 3. Top: V = 0.5. 
Bottom: = 1.0. 

charge fluctuations play only a minor role compared to 
the SDW phase, similar to the one-dimensional system. 

VI. CONCLUSIONS 

In this paper we have presented a generalization of 
the variational cluster perturbation theory to extended 
Hubbard models at half filling. The method is based on 
the self-energy-functional approach (SFA) which uses dy- 
namical information of an exactly solvable system (ref- 
erence system H') in order to approximate the physics 
in the thermodynamic limit. For the application of this 
method, a mean-field decoupling of the inter-cluster part 
of the nearest-neighbor Coulomb interaction is performed 
first. After this step, one is left with a Hamiltonian which 
couples the different clusters via the hopping only and 
which can be treated by the known (variational) CPT 
procedure. The mean-field decoupling yields effective 
onsite potentials on the cluster boundaries as external 
parameters of the Hamiltonian. These parameters are 
determined either self-consistently on an isolated cluster 
(sufficient for the study of first order phase transitions) 
or by determination of the minimum of the SFA grand 



potential. 

In order to test the accuracy of our approach we ap- 
plied the method to the extended Hubbard model in one 
dimension, because results from other methods like QMC 
and DMRG are available for comparison. At J7 = 8 
the results for the critical interaction Vc, the ground- 
state energy, kinetic energy, and charge order param- 
eter showed excellent quantitative agreement with pre- 
vious QMC studies. At ?7 = 3 our method predicted 
a second-order phase transition with transition point 
Vc — 1.665(5) again in good agreement with previous 
studies. The ground state energy, kinetic energy, and 
order parameter do not seem to have been calculated be- 
fore. 

In addition we calculated the spectral function for sev- 
eral values of the interaction V , which has not been done 
previously. At both U — % and U = 3, wc found evi- 
dence for spin-charge separation in the SDW phase, but 
not in the CDW phase. By fitting the bands by Hartree- 
Fock dispersions we found that the hopping parameter 
is strongly renormalized. The agreement between the 
fitted value of the gap and the value within the Hartree- 
Fock approximation was much better in the CDW phase 
than in the SDW phase giving rise to the conclusion that 
charge fluctuations play a minor role in the CDW phase. 

Whereas the application of sophisticated methods like 
DMRG or fcrmionic loop-update QMC to more than 
one dimension is difficult, this extension is straightfor- 
ward within the present approach. We were thus able 
to perform the first non-perturbative study of the two- 
dimensional extended Hubbard model on a square lat- 
tice at half filling and zero temperature beyond Hartree- 
Fock. We found first order transitions at both [/ = 8 
and [/ = 3 with transition points Vc — 2.023(1) and 
Vc — 0.770(3) for an Nc = 48 super cluster, respectively. 
The spectral function in the SDW phase shows coherent 
low-energy quasi-particle excitations with band width set 
by the magnetic exchange constant J, and an incoherent 
background, consistent with previous QMC studies for 
the Hubbard model at F = 0. The Hartree-Fock pre- 
diction differs significantly from the low-energy feature 
and does not describe the splitting into coherent quasi- 
particle bands and incoherent background. In the CDW 
phase the Hartree-Fock dispersions account much better 
for the excitations, and no additional low-energy features 
caused by a magnetic origin could be found. Similar to 
one dimension the agreement between the Hartree-Fock 
approximation and the low-energy excitations obtained 
by the present method is much better in the CDW phase, 
confirming that charge fiuctuations are less important in 
the charge-ordered phase than in the SDW phase. 

In this paper we applied our method to half-filled sys- 
tems only, but one can study ordering phenomena at 
other fillings, too, as long as the possible order patterns 
are commensurate with the shape of clusters used as ref- 
erence system. With some effort it is also possible to 
study phases with long wave-length charge density waves 
by coupling several clusters to a super cluster and ap- 
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plying appropriate continuity conditions between the in- 
dividual clusters within the super cluster. In addition 
the application to systems with lattice geometry differ- 
ent from the two-dimensional square lattice, e.g. ladder 
materials, is an interesting subject for further studies. 
Work in this direction is in progress. 
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where we assumed without loss of generality that the 
bonds [ij] connect sites i on sublattice A with sites j on 
sublattice B. The free energy of the system is given by 



F = - 



4lnZ = -ilntre-'^^"''^'') 



- — In tr exp 



13' 



R 



(A3) 



Taking the derivative with respect to S yields 



nRi - riRj + 26] 




(A4) 



APPENDIX A: 



MEAN-FIELD SOLUTION AND 
FREE ENERGY 



In this section we show that a mean-field solution ob- 
tained self consistently is directly connected to a min- 
imum in the free energy. For simplicity let us assume 
that we have only two different mean-field parameters 
Ayi = 1 — (5 and Xb = I + S, see also Sec. lIIII We can 
write the mean-field decoupled Hamiltonian Eq. 11()|) as 



E 

R 



(Al) 



where includes all terms independent of the mean- 
field parameters. According to the third line in Eq. (|13|1 . 

■f^MkR-'*^) is given by 

M 

^VJ2 + + S^)] (A2) 



All clusters are equivalent, therefore we suppress the in- 
dex R in the following. Setting this derivative to zero we 
get the self-consistency condition 



m 

For one dimension, Eq. (|A5|I is given by 

{tin) - (ni) = 2(5, 



(A5) 



(A6) 



because in this case we have only one decoupled bond 
[liV] with site 1(N) belonging to sublattice A{B), re- 
spectively. To conclude, one can state that if self con- 
sistency. Eg. (jA5|) . is fulfilled, then the free energy has 
an extremum with respect to the mean-field parameter 
d. By thermodynamic stability arguments this extremum 
always has to be a minimum. 
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